/*  sdcchpi.cpp  by K.Tsuru  */
// ver. 2.30
// function ID 3555   DRADIX
/***************************************************************
SN library
An evaluation of pi by the binary splitting method using Chudnovskys' formura.

 1             1                           \prod_{k=0}^(n-1) {-(2k+1)(6k+1)(6k+5)}
---- = -------------------sum_{n=0]^{L}-------------------------------------------- (A+B*n)
 pi    426880*sqrt(10005)                  \prod_{k=0}^n {(C*k)^3}/24 (k>=1)

where  A = 13591409, B = 545140134, C = 640320 and C^(1.5)/12=426880*sqrt(10005). See below about "L".
Let
 426880*sqrt(10005)
--------------------- = sum_{k=0]^{L}P(k)
        pi

P(k+1)   -24(6k+1)(2k+1)(6k+5)(A+Bk+B)
------ = ---------------------------
 P(k)      C^3(k+1)^3(A+Bk)
                   P(k+1)
\lim_{k \to \infty}------ = -24*(6^2*2)/C^3 = -1/15193 13730 56000 = -1/Q
                    P(k)
(i)log10(Q) = 14.18...

P(1)                                A+B
---- = -120*(A+B)/(A*C^3) = - ------------------- = - 1/Q0
P(0)                           2187811772006400*A
                           13591409
(ii)Q0 = 2187811772006400----------- = 5321 95559 40398.61289....
                          558731543
    log10(Q0) = 13.726...
Then adding one term the precision raises up about 14 digits.
The upper limit "L" may be chosen L= (the effective number of digits)/14.

---- 32bit SLong version ----------------
Binary Splitting method Nov.19,2016
CPU Core i7-3770 3.4GHz Memory 32.0GB
 digits   | CPU time(sec)
  4000069 |  21.0
  8000029 |  46.75
 33550029 | 233.121 (33 million)
 50000029 | 372.56[6min 12sec]
100000029 |1140[19min]+output= 12(sec)(100 millon)
---- 64bit SDouble version ----------------
 digits   | CPU time(sec)
 500000   | 1.78949(sec)
 50000028 | 442.53(sec) [7.3755 min]
***************************************************************/

#ifndef SN_H
#include "sn.h"
#endif

static const double upPerTerm = 14.18;
#define CH_PI_SLONG 0
#if CH_PI_SLONG
////// SLong version ///////  1.75(sec) 500029 digits

class BSChudnovskysPi : public BinarySplitting<SLong> {
public:
    /* Constructer : 'L' is upToTerm, 'f' is precision*/
    BSChudnovskysPi(long L, long f) : BinarySplitting<SLong>(L, f){}

    void setABC(long k, SLong& a, SLong& b, SLong& c) {
        static const SLong tA(13591409L);
        static const SLong tB(545140134L);
        static const SLong tC(640320L);
        static const SLong tC3((tC * tC * tC)/24L);

        a = tA + tB * k;
     // B(k) = - (2 * k + 1) *(6 * k + 1) * (6 * k + 5)
        long p = - (2L * k + 1L), q = 6L * k + 1L, r = 6L * k + 5L;
        b = p;  b *= q;  b *= r;
    // C(0) = 1, C(k) = tC3 * k^3(k >= 1)
        if(k == 0) c = 1;
        else {
            c = k;  c *= (c*c);  c *= tC3;
        }
    }

    SDouble getValue() {
        putTogether();
        SDouble pi = BinarySplitting<SLong>::getValue(true);
        A=0.0;  B= 0.0; C=0.0;  // free memory

        SDouble a = Sqrt(10005L);
        a *= 426880L;
        pi *= a;
        return pi;
    }
};
SDouble ChudnovskysPi() {
    SDouble tmp;
    long f = long(tmp.EffFig() + tmp.Hidden()+ 1u)*DFIGURES;
    long L = (long)((double)f/upPerTerm);

    BSChudnovskysPi pi(L, f);
    return pi.getValue();
}

#else
////// SDouble version /////// 4.92(sec) 500029 digits

class BSChudnovskysPi : public BinarySplitting<SDouble> {
public:
    BSChudnovskysPi(long L, long f) : BinarySplitting<SDouble>(L, f){}

    void setABC(long k, SDouble& a, SDouble& b, SDouble& c) {
        static const SDouble A(13591409L);
        static const SDouble B(545140134L);
        static const SDouble C(640320L);
        static const SDouble C3((C * C * C)/24L);
        a = A + B * k;
     // B(k) = - (2 * k + 1) *(6 * k + 1) * (6 * k + 5)
        long p = - (2L * k + 1L), q = 6L * k + 1L, r = 6L * k + 5L;
        b = p;  b *= q;  b *= r;
    // C(0) = 1, C(k) = C3 * k^3
        if(k == 0) c = 1;
        else {
            c = k;  c *= (c*c);  c *= C3;
        }
    }

    SDouble getValue() {
        putTogether();
        SDouble pi = BinarySplitting<SDouble>::getValue(true), a = Sqrt(10005L);
        pi *= a;
        pi *= 426880L;
        return pi;
    }
};

SDouble ChudnovskysPi() {
    SDouble tmp;
    long f = long( tmp.EffFig() + tmp.Hidden() ), L = (long)((double)f / 2);

    BSChudnovskysPi pi(L, f);

    return pi.getValue();
}
#endif